#ifndef MAT4_H_
#define MAT4_H_

#include "../common.h"

#include "Mat3.h"

namespace physicsworks {

namespace math {

template<class T, class Tr = real_traits<T> > class Mat4 : public Mat<T, Tr> {
private:
	typedef Mat<T, Tr> Base;

	/**
	 * Returns a 3x1 matrix as Vec3<T> supplying the column number.
	 * This matrix is formed by Aic for 0 <= i < 3.
	 */
	Vec3<T> vector3(size_t c) {
		Vec3<T> v;
		for (size_t i = 0; i < 3; ++i)
			v(i) = Base::operator()(i, c);
		return v;
	}
	
public:

	using Mat<T, Tr>::operator=;

	Mat4<T, Tr>() : Mat<T, Tr>(4, 4) {}
	
	Vec3<T> x() const { return vector3(0); }
	Vec3<T> y() const { return vector3(1); }
	Vec3<T> z() const { return vector3(2); }
	Vec3<T> w() const { return vector3(3); }

	T determinant() const;
	
	void genTransformation(const Vec3<T>& position, const Mat3<T>& orientation) {
		std::valarray<T>& arr = Base::arr;
		arr[0]  = orientation(0,0);	arr[1]  = orientation(0,1); arr[2]  = orientation(0,2); arr[3]  = position.x;
		arr[4]  = orientation(1,0); arr[5]  = orientation(1,1); arr[6]  = orientation(1,2); arr[7]  = position.y;
		arr[8]  = orientation(2,0); arr[9]  = orientation(2,1); arr[10] = orientation(2,2); arr[11] = position.z;
		arr[12] =                0; arr[13] =                0; arr[14] =                0; arr[15] =          1;
	}
};

/**
 * Returns the matrix's determinant.
 */
template<class T, class Tr> T Mat4<T, Tr>::determinant() const
{
	Mat4<T, Tr>& m = *this;
	
	return  ( m(3,3)*m(4,4) - m(3,4)*m(4,3) ) * ( m(1,1)*m(2,2) - m(1,2)*m(2,1) )
	       +( m(3,4)*m(4,2) - m(4,4)*m(3,2) ) * ( m(1,1)*m(2,3) - m(1,3)*m(2,1) )
	       +( m(3,4)*m(4,1) - m(4,4)*m(3,1) ) * ( m(1,3)*m(2,2) - m(1,2)*m(2,3) )
	       +( m(4,3)*m(3,2) - m(3,3)*m(4,2) ) * ( m(1,1)*m(2,4) - m(1,4)*m(2,1) )
	       +( m(3,3)*m(4,1) - m(4,3)*m(3,1) ) * ( m(1,2)*m(2,4) - m(1,4)*m(2,2) )
	       +( m(4,2)*m(3,1) - m(3,2)*m(4,1) ) * ( m(1,3)*m(2,4) - m(1,4)*m(2,3) );
}

}

}

#endif /*MAT4_H_*/
